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PREFACE 


Satellite-based altimetric data taken by GEOS-3 and SEASAT over the Black Sea and Caspian 
Sea are analyzed and a least-squares collocation technique Is used to predict the geoid 
undulations on a .25° x .25® grid and to transform these geoid undulations to free air gravity 
anomalies. Rapp's 180 x 180 geopotential model is used as the reference surface for the 
collocation procedure. The result of geoid to gravity transformation is. however, sensitive to 
geographic variability in the information content of the reference geopotential model used. 

For example, detailed surface gravity data are available for incorporation into the reference 
model over the Black Sea, resulting in a reference model with significant Information content 
at short wavelengths. It is shown that reliable estimation of gravity anomalies from gridded 
geoid heights is generally possible over regions such as the Black Sea. using the conventional 
collocation technique. Over regions where surface data are generally not available for 
incorporation into the reference model, such as the Caspian Sea, an enhanced algorithm is 
needed to obtain reliable gravity anomalies. A feasible algorithm would extract both short and 
long wavelength information from the altimetric geoid heights and create an accurate gravity 
field at all wavelengths. Methods of generating such an algorithm are described and tested. 
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L INTRODUCTION 


This project entailed processing satellite altimeter data over inland seas (the Black and 
Caspian Seas) for recovery of area mean gravity information. Gravity Information in this area 
of the world is not readily available, so the possibility of obtaining it from the processing of 
altimeter observations is attractive. The mean surface level of the seas approximates an 
equipotential surface. Therefore, information about the underlying gravity potential and its 
derivative, gravity, can be obtained from measurement of the relative shape of this surface by 
means of altimetry. 

Our principal objective on this project was to complete and extend analyses done in a 
previous study, verify those results, and document the results and techniques. A secondary 
objective was to improve the algorithms and results, if possible. The basic approach used by 
the previous study, and followed by STX personnel, involved five steps: 

1. Edit geoid height data to remove overland data; 

2. Evaluate geoid height differences at crossover points; 

3. Remove orbit errors from geoid heights using crossover differences; 

4. Grid geoid height data at 0.25° x 0.25° Intervals; 

5. Estimate 0.25° x 0.25° gravity anomalies from gridded geoid heights using the 
collocation technique. 

The need for step 1 is obvious. Steps 2 and 3 are necessary because satellite altimeter 
measurements cannot yield accurate sea surface heights unless differences in satellite heights 
(orbit differences) from pass to pass are rectified and reduced to a common reference. If the 
mean sea surface height at a given location is constant over the time span of the altimeter 
data used, any difference in surface height between two altimeter passes at the point where 
they cross Is due to orbit differences (differences of up to about 50 cm could be a consequence 
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of tides, especially solid earth tides, whose amplitudes are about 25 cm). Since the orbit 
differences are constant for the short arcs over the inland seas, removal of a constant bias 
from all the crossover differences of a given pass should effectively rectify orbit differences. 
Area mean surface height values are determined and reduced to the reference geoid in step 4. 
In step 5 these area mean geoid heights are processed and area mean gravity anomaly values 
are predicted using a linear least-squares estimation technique, called collocation, formulated 
by Moritz ( 1978). The collocation technique Is essentially a differential operation transforming 
geopotential information to its first derivative, gravity. Knowledge of the statistical correlation 
between area mean geoid heights and gravity anomalies Is required in the geoid-to-gravity 
transformation. 


GEOS-3 altimeter data is of lesser quality (standard deviation between 25 and 50 cm, 
depending on operating mode) than that of SEASAT (6 to 10 cm), primarily because SEASAT 
uses an advanced radar altimeter design. The GEOS-3 altimeter operating modes, intensive and 
global, are differentiated primarily by data rate, which explains the corresponding difference 
in quality. The GEOS-3 mission collected data between 1975 and 1978 over latitudes up to 65 

degrees, whereas SEASAT collected data only during 100 days in 1978 over latitudes up to 72 
degrees. 


In the next section, steps 1 to 4 are discussed in greater detail. An expanded 
mathematical description of the collocation technique and results of its application to both geoid 
gridding and gravity prediction are presented in Section III. The results of this project are 
discussed in Section IV. Conclusions and recommendations are presented in Section V. 
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IL ALTIMETER DATA 


The altimeter data over the Black and Caspian Seas, obtained from NASA / GSFC, had been 
processed with GEODYN and written in GEODYN format Program PRTDATA (Au et at., 1989) 
organized the data records and extracted relevant geodetic and auxiliary Information. There 
are 83 GEOS-3 and 62 SEASAT passes over the Black Sea written in 9,248 data records. Over 
the Caspian Sea there are 71 GEOS-3 and 23 SEASAT passes written in 20,642 data records. 
Progr am PLTGRP (Au et at., 1989) was developed to plot the altimetric surface elevation of 
each pass as a function of time. Typical elevation profiles over the Black and Caspian Seas 
are shown In Figures 1 and 2. respectively. Visual examination of these surface elevation 
profiles suggests that the data are relatively noiseless, except for a few occurrences of data 
spikes and data gaps. In subsequent data processing programs, data spikes were eliminated by 
removing data points that deviated from adjacent ones by more than a given value, which was 
2 meters for the Black Sea and 5 meters for the Caspian Sea. The ground tracks of these 
satellite passes over the Black and Caspian Seas, plotted with Program GRNTRK (Au et at., 
1989), are shown in Figures 3 and 4, respectively. 


Crossove r Adjustments 

The major error source in altimetric geoid undulations is the uncertainty in the radial 
component of the satellite trajectory. This uncertainty is manifested in the misclosure of 
surface elevation at ground track intersections (crossovers) between passes. For the short arcs 
considered here, the orbit error can be modeled as a bias applied to all the data of a given 
pass. The optimum biases are such that crossover differences are minimized, holding one pass 
fixed so that all the satellite passes can be defined with respect to a common reference model. 
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Figure lb. Sample altimeter elevation profiles of the Black Sea. 
Profile with data spike. 
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Figure 2a. Sample altimeter elevation profiles of the Caspian Sea. 
Relatively noiseless profile. 
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Figure 2b. Sample altimeter elevation profiles of the Caspian Sea. 
Profile with data gap. 
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Figure 3. GEOS-3/SEASAT ground tracks over the Black Sea. 
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Figure 4. GEOS-3/SEASAT ground tracks over the Caspian Sea. 
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To calculate crossover differences, one must first locate the crossover location in latitude 
and longitude. There are various methods by which this point can be determined. We have 
adopted an analytical method of modeling the ground track of the pass. For relatively short 
arcs, such as the satellite passes over the Black and Caspian Seas, the ground track can be 
approximated by a second-degree equation, 

Y « aX 2 + bX + c j 

where Y and X are. respectively, the latitude and longitude vectors of ground track records, 
and a, b and c are polynomial coefficients to be determined by fitting the ground track data 
using the method of least-squares. The error in this satellite arc representation is less than 
the radius of the illuminated area of the altimeter signal at the sea surface. The latitudes of 
the crossover point of two passes, Y ! and Y 2 , are set equal. V , 1 = Y 2 . The longitude at which 
this crossover occurs is determined by solving the quadratic equation for X. An acceptable 
crossover point, naturally, must lie within the ground track records of both passes considered. 
Program XOVERO (Au et aL, 1989) was used to determine the latitudes and longitudes of all 
possible crossovers and to prepare for subsequent crossover adjustments. 


Once the crossover point is located, the altimetric height Is Interpolated by cubic splines 
from the nearest data for each pass. The true geoid undulation at a crossover point must be 
the same for both passes regardless of satellite and time. Altimetric height, however, is not 
exactly the same as geoid undulation. For example, temporal processes such as solid earth and 
ocean tides may cause the sea surface height to be different at the different times of the 
crossing altimeter passes. Ocean tides on small seas like these should contribute less than 10 
cm to the crossover difference, but diurnal earth tides may be expected to contribute up to 
about 80 cm. Fortunately, tides are such broad-scale features in both space and time that they 
are manifest as a constant bias in a single pass of altimeter measurements over features as 
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small as the Black and Caspian Seas. Thus any earth tide effects alias with the orbit error 
bias and are removed when this bias is adjusted. To Illustrate this adjustment, let H ( ° be the 
true geo id undulation for pass ( at a crossover point and b ( be the bias assumed for this pass. 
The observed geoid undulation Hj is given by 

H t - H ( ° + b, . (2) 

The difference dy at a crossover point between pass l and pass J will be 

d t j = Ht-Hj 

= ( H { ° + b t ) - { Hf + bj ) 

= b t -bj (3) 

because H t ° and Hj° must be identical at a crossover point. An overdetermined system of 
equations in b results if all crossover residuals over an inland sea are considered. The bias 
for each pass, therefore, can be determined using the method of weighted least-squares, thus 
minimizing the crossover differences d. The standard error is assumed to be 25 cm for GEOS- 
3 and 10 cm for SEASAT in the weight matrix. The optimal pass bias 
vector B is given by 

B = {AWAr l lAWD) (4) 

where D is the vector of crossover differences, W is a diagonal matrix in which diagonal 
elements are the sum of the inverse of the squares of the assumed sigma for each satellite. 

The matrix A is sparse. Each row of A contains all zeros except unitary value in the column 
associated with a pass t and a negative unitary value in the column associated with pass J . 
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The pass with the most crossovers was chosen as the reference pass. The bias for this 
pass Is not estimated, but, after the crossover adjustment process, is assigned the average geoid 
height along the reference pass ground track, as calculated from the reference geopotential 
model. Rapp’s 180 x 180 model is the reference geopotential model used in the current report 
Program XOVER (Au et aL, 1989) was developed to perform the crossover adjustment process. 
An error covariance matrix of the crossover adjustment was also determined by XOVER. This 
error covariance matrix adds to the error associated with each satellite pass a crossover 
adjustment error. 

From the geometry of altimeter passes in the current dataset, there are at most, 2109 
crossovers over the Black Sea, 570 of which are GEOS-3 with GEOS-3, 383 are SEASAT with 
SEASAT and 1 156 are GEOS-3 with SEASAT. These possible crossover locations were carefully 
checked in order to eliminate any crossover locations from further consideration If they 
occurred at data gaps, which are defined to be part of a satellite arc that did not have an 
altimeter observation within 70 km (about 10 seconds in time). Such editing reduced the 
number of crossovers to 1891, 521 of which are GEOS-3 with GEOS-3, 350 are SEASAT with 
SEASAT and 1020 are GEOS-3 with SEASAT. The RMS (root-mean-square) of the crossover 
residuals before crossover adjustment is 3.96 m. The RMS after crossover adjustment is 
reduced to 25 cm. The reference pass is the GEOS-3 pass #10557. Typical adjusted profiles 
and their corresponding reference models are shown in Figure 5. Program PLTAEP (Au et ai f 
1989) was developed to plot the adjusted profiles. 

We have analyzed these crossover difference statistics for consistency with the reported 
precision of the GEOS-3 and SEASAT altimeter data. Stanley (1979) reported the standard 
deviation of the GEOS-3 altimeter precision at 25 cm in the intensive mode, and 50 cm in the 
global mode. We have no knowledge of the mode of the GEOS-3 altimeter data used in this 
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Figure 5a. Sample adjusted altimeter elevation profiles of the Black Sea. Rapp’s 180 x 180 
geoid model is represented by the dash curve. 

Profile of the reference pass. 
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Figure 5b. Sample adjusted altimeter elevation profiles of the Black Sea. Rapp s 180 x 180 

geoid model is represented by the dash curve. 

Typical profile. 
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study, because It has been reprocessed Into a consistent format that obscures the mode. 
However. Wagner (1979) found a total RMS of 33 cm in the high-frequency spectra from 47 
passes of global and intensive GEOS-3 altimeter data. For the SEASAT data, Townsend (1980) 
reported a data noise level between 6 and 10 cm, with the higher values occurring at higher 
sea states. 

The RMS of crossover differences consists of the square root of the mean of a weighted 
sum of squares of data of different precisions. If we assume that the data of the two crossing 
passes is uncorrelated, then the sum of squares is comprised of three terms: 

1) twice the variance of GEOS-3 data, weighted by the number of GEOS-3 crossovers: 

2) twice the variance of SEASAT data, weighted by the number of SEASAT crossovers; 

3) the sum of variances of GEOS-3 and SEASAT data, weighted by the number of mixed 
crossovers. 


If we insert the above satellite data variances into the above sum-of-squares equation, we may 
establish the reality of our observed crossover RMS. Inserting [25 cm| 2 and [ 10 cmP for 
GEOS-3 and SEASAT variances, respectively, an RMS crossover difference of 0.28 m is predicted, 
which is in good agreement with our crossover RMS of 0.25 meters. 

Over the Caspian Sea there are 1074 possible crossovers. 539 of which are GEOS-3 with 
GEOS-3, 79 are SEASAT with SEASAT and 456 are GEOS-3 with SEASAT. The crossover- 
cleanup process described above for the Black Sea was also applied to the Caspian Sea 
crossovers. Two passes that cross no other passes were removed from the dataset. This 
reduced the number of crossovers to 972, of which 494 are GEOS-3 with GEOS-3, 77 are 
SEASAT with SEASAT and 401 are GEOS-3 with SEASAT. The RMS of crossover differences 
before crossover adjustment is 2.21 m. The RMS after crossover adjustment is reduced to 37 
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cm. The reference pass is the SEASAT pass #832. Typical adjusted profiles and their 
corresponding reference models are shown in Figure 6. 

Following the method described above, if we insert reported variances for data of one 
satellite into the expression for the RSS (residual-squared sum) of the Caspian Sea crossover 
differences, the resulting variance of the other satellite’s data may be checked against reported 
precision values. If, for example, we insert a variance corresponding to 6 cm noise for SEASAT, 
we find that the corresponding GEOS-3 noise level satisfying the crossover RMS is 32 cm. If 
we Insert a variance for SEASAT corresponding to a 10-cm noise level, the corresponding GEOS- 
3 noise level satisfying the crossover RMS is about 31 cm. Both of these numbers are quite 
reasonable according to the references cited above. On the other hand, If we insert a GEOS- 
3 noise level of 25 cm, we find that this implies a SEASAT noise level of 31 cm, which is 
unreasonably high, whereas inserting a GEOS-3 noise level of 50 cm is prohibited (requires 
square roots of negative numbers). Thus our crossover RMS of 0.38 meters is consistent with 
reasonable values of SEASAT and GEOS-3 altimeter precisions. 

After the data were corrected for pass biases, an overall bias representing the average 
difference in height between the reference pass and the reference geoid is added to the data. 

For the Caspian Sea, the adjustment was -34 m, and for the Black Sea data, 0 m. 


Preliminary Analysis of t he Adjusted Data 

The adjusted inland sea. data must be gridded for geodetic collocation analyses. This is a 
time-consuming process unless the data records are properly arranged. To group data records 
by proximity to the same grid point, program SORT (Au et al., 1989) reads through a dataset 
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Figure 6a. Sample adjusted altimeter elevation profiles of the Caspian Sea. Rapp s 180 x 180 

geoid model Is represented by the dash curve. 

Profile of the reference pass. 
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Figure 6b. Sample adjusted altimeter elevation profiles of the Caspian Sea. Rapp s 180 x 180 

geoid model Is represented by the dash curve. 

Typical profile. 
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and determines the grid points, depending on the chosen cap size, with which a particular data 
record is associated. Identification indices are then assigned to each data record so that data 
records associated with a common grid point can be grouped together by a sort / merge process. 
In the current report, the grid size is a quarter degree. The cap radius chosen for subsequent 
data analyses is also a quarter degree, so that all data within a quarter degree of a given grid 
point are associated with that grid point. Note that it is possible for a given data record to 
be associated with as many as four grid points. 

Geoid undulations of the Inland seas were initially gridded according to the method of 
weighted averages using Program WGTAVG (Au et al„ 1989). The weight of a data point with 
respect to a grid point is a function of the square of its distance from the grid point. Contour 
maps of the weighted-average geoid undulations of the Black and Caspian Seas are shown in 
Figures 7 and 8. respectively. The contour maps of the geoid undulations of the Black Sea 
and Caspian Sea according to the reference geopotential model, Rapp’s 180 x 180 model, are 
shown in Figures 9 and 10, respectively. Program CONTOUR (Au et al., 1989) is used for 
contouring both the gridded data and reference models. 

The contour maps of the weighted-average and reference geoid undulations of the Black 
Sea closely resemble each other. The geoid contour maps of the Caspian Sea. however, show 
substantial differences. Assuming comparable short-wavelength information content in the 
altimeter data over both seas, it can be concluded that there is more short-wavelength 
information in the reference geoid undulation model of the Black Sea than in the model of the 
Caspian Sea. The information content of a reference geopotential model has considerable effect 
on the quality of gravity anomaly results transformed from geoid undulation data, as will be 
shown in later sections. 
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Figure 7. A contour map of the weighted-average geoid undulations 
(m above mean sea level) of the Black Sea. 


20 





Oft/Gffl.VX p* » 

° F PQQ'i) p //," 



Figure 8. A contour map of the weighted-average geoid undulations 
(m above mean sea level) of the Caspian Sea. 
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Figure 10. A contour map of Rapp's 180 x 180 reference geoid undulations 
(m above mean sea level) of the Caspian Sea. 
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IIL BRIEF DESCRIPTION OF LINEAR LEAST-SQUARES INTERPOLATION AND COLLOCATION 


Collocation is a predictive method based on linear least-squares interpolation. A general 
description of collocation can be found in Moritz ( 1978). A brief review of the collocation 
method Is given in this section. 

LetTTP) be a predicted value of a function at a point P. The predicted value can be expressed 

by a linear combination ofm observable (stochastic variables) in the neighborhood of P, such 
that 


m 

J\P) X ^pi^pi 
<=1 

- XpNp 


where N p is an m * 1 vector of observables and a p is a set of nonunique coefficients. 
Is the true value of the function at P. the local interpolation error * p at P is given by 


(5) 
if r p 


fp = T p - (TpNp 
= { 1 -ffp } 



( 6 ) 


Similarly, interpolation error at another point Q is 


e Q ~ t * -&Q ) 
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The error covariance for points P and Q is 
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Let M be an average operator such that 
opg = M { «pg) 

Is the error covariance of prediction errors <p and < g at points P and 9- A covariance 
matrix, I , can be defined as 


1 = M 


T'P 

Np 


{ T Q "9 > 


>9 

C PN 

~N9 

C NN 


( 10 ) 


where C pg is the error covariance of points P and 9 . Cptf *» an n x 1 covariance vector 
relating the elements of Ng to Tp, C wg Is an m * 1 covariance vector relating the elements 
of Np to T g . and C N fi is an m * n covariance matrix of the elements of N P and Ng. 

Elements of the covariance matrix are values of a chosen covariance function that reflect the 
interdependence of two observables. A covariance function, therefore, can be any mathematical 
function, e.g., an ensemble average of observables, or an empirical formula. 


The error covariance of P and Q is then given by 


(Jpg = Cpg - XgCp N - <rpC QN + ^pC^/dp 
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and the variance at P is 


°P L C o ~ 2ff P c PiV + a P c NN a P 
where C Q is the assumed error at the prediction point P. 


( 12 ) 


The necessary condition for a minimum of of is that 
2 

top / dcip — 0 

r y ( 13 ) 

which implies that 

°p = {C nnT 1c pn ■ 

Substituting a p into the equaUons for J{P) and for the variance at P, the predicted value 
at P and Its variance are given by 

JiP) = c PN (c N fjr l N p (i5) 

= c pp - CpNiCflfir l c PN ( 

The covariance matrix previously derived is only for the stochastic variables. If there are 
random errors associated with the stochastic variables, the stochastic covariance matrix should 
be supplemented by an error covariance matrix D N fj, such that 

C (( = C NN + °NN (17) 

Is the covariance matrix. The predicted value at P becomes 

JXP) = C pN C (e ~ l Np (lg) 

where C pN remains unchanged. The notation that indicates a transpose. ~ . will be omitted in 
the following sections. 
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IV. APPLICATION OF COLLOCATION TECHNIQUE 


r.ni location O Hdding of Geoid Undulatlgoa 

According to the linear least-squares Interpolation formula, the predicted geold undulation 
NIP) at a point P Is given by 

N(P) = C pN {C NN r l Np 

where C pN is a covariance vector relating the undulation at P to the observable in the 
neighborhood of P. and C NN is the covariance matrix. The covariance matrix is the sum of 
the stochastic covariance matrix and the error covariance matrix. In the current report, the 
error covariance matrix Is a combination of the random error associated with each observable 
and the error covariance matrix obtained from crossover adjustments. The stochastic 
covariance matrix is derived from a local residual covariance function based on the difference 
between gridded weighted-average geoid undulations and the reference geopotential model. A 
convolution technique described by Moritz ( 1978) is used to determine the local covariance 
function. The local covariance function is constructed by convolution of the difference 
between the weighted-average geoid data and the reference geoid. The resultant covariance 
function is, in effect, a least-squares niter (Treitel and Robinson, 1966). This filter determines 
the contribution of each observable to the predicted value at a grid point. It is observed that 
the interpolation results are rather insensitive to the covariance function being used as long 
as reasonable effective weights are assigned to each observable. For example, if the weights 
are assigned as a function of the inverse of the square of distance from the point at which 
prediction is made, the weighted-average results coincide with the collocation results. 
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The variance at each gridded point is given by 

a 2 =c Q - c PN ic NN r 1 c PN 


( 20 ) 


where C Q is the assumed value of the covariance function at the predicted point. 


Program GEOID (Au et al., 1989) was used to perform the geoid undulation Interpolation, 
and Program EMPCOV (Au et al., 1989) was used to determine local empirical covariance 
functions. The input dataset used in Program GEOID is a sorted version of the output from 
Program SORT, in which cutoff cap size is set to one quarter degree. The effective integration 
cap size of this interpolation is, therefore, also a quarter degree. The Integration cap radius 
effectively excludes longer wavelength information from the gridding. 

Contour maps of the geoid undulations of the Black and Caspian Seas gridded by the 
collocation method are shown, respectively, In Figures 1 1 and 12. The square root of the 
variance of the gridded data of the Black and Caspian Seas are plotted, respectively, in Figures 
13 and 14. Note that this error in the gridded data varies Inversely with the density of ground 
tracks (Figures 3 and 4) and is much smaller ( 10-cm RMS) for the Black Sea than for the 
Caspian Sea (20-cm RMS). This disparity in error values reflects somewhat the larger crossover 
error In the Caspian, but mostly the less uniform and sparser data distribution. Plots of the 
Black Sea and Caspian Sea local empirical covariance function are shown in Figures 15 and 16, 
respectively. Note that the weighted-average and collocation geoid undulation maps are very 

similar, demonstrating that collocation, when used as an interpolator, is not sensitive to the 
choice of covariance function. 
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Figure 1 1. A contour map of the collocation geoid undulations 
(m above mean sea level) of the Black Sea. 


29 







ORIGINAL PAGE IS 

OF POOR QUALITY 


Figure 12. A contour map of the collocation geoid undulations 
<m above mean sea level) of the Caspian Sea. 
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Figure 13. A contour map of the square root of the variance (cm) of the 
collocation geoid undulations of the Black Sea. 
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Figure 14. A contour map of the square root of the variance (cm) of the 
collocation geoid undulations of the Caspian Sea. 
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Figure 15. Normalized local residual empirical covariance functions of the Black Sea. 
The normalization coefficients for: 

1) undulation-undulation (solid line) is 1.61 m 2 , 

2) undulation-gravity (dash line) is 10.00 m-mgal, and 

3) gravity-gravity (dotted line) is 125.5 (mgal) 2 . 
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ngurc 16. Normalized local rea.dual empirical covariance redone or ,he Caeplan Sea 
The normalization coefficients for: 

1) undulation-undulation (solid line) is 8.48 m z, 

2) undulation-gravity (dash line) is 42. 15 m-mgal. and 

3) gravity-gravity (dotted line) is 340.2 (mgalp. 
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F.fltimat i pn of Gravity Anoma li e s 


Following the collocation method deecrlhed above, geold undulation, can be transformed 
into gravity anomalies (Rapp. 1986) according to the equation 

i * (21) 

Ag = Cg^iC^+DF IN-N R ) + &g R 

where Ar, I. the predimed point gravity anomaly. C gN 1. the covariance vector of geold-to- 
gra.lt, traneformation, C m le the covariance matrix for undulatlon-undulatlon interpolation. D 
the error covariance matrix that la conatructed In part from the variance of the prevloua 
Interpolation of geold undulation., N la the vector of gridded geold undulation.. Nr la the 
vector of reference model geold undulation, that correaponda to each obaerved value of N. and 
4 9r „ the reference model gravity anomaly value at die predicted grid point. The variance .a 

given by 

. ( 22 ) 
<7* = Cgg - C gN {C NN +Dr C gN 


where C gg la the aaaumed value of the covariance function at the predicted grid point. 

Program GRAVAN ,Au « nl.. 19891 wa. uaed to perform the geold-to-gravlty traneformation. 

The geold undulation., gridded by collocation technique, and the correapondlng variance., were 
uaed as Input dataaeta. The demand for computer resource, was greatly reduced by the choice 
of gridded data because a gridded dataset contain, far fewer data point, than an ungrldded 

dataset. 


In the current study, three sets of covariance functions were investigated for the 
undulation-to-gravity transformation. They are: 1, Rapp s empirics, covariance functions based 
on a 180 x 180 global geopotential model, 2) theoretical self-consistent covariance functions 
(Jordan, 1972), and 3) local residual empirical covariance functions derived from iterative use 
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of the convolution method previously described. In some cases, it was necessary to augment 
the a-prlort error matrix D In order to achieve stability in the gravity prediction. For the 
local empirical covariance functions. [60 cm)* of variance was added to the diagonal of D. 
whereas for the Rapp and Jordan models, only [25 cm) 2 of variance was added. 

Plots of Rapp’s 180 x 180 global covariance functions are shown In Figure 17. Contour 
maps of the estimated Black Sea gravity anomalies and the square root of their variance, using 
Rapp s covariance functions, are shown, respectively. In Figures 18 and 19. Corresponding 
results for the Caspian Sea. again using Rapp’s covariance functions, are shown in Figures 20 
and 21. As with the gridded geoid undulation error, the pattern of gravity anomaly error 
correlates inversely to the data distribution, but the magnitude of error is nearly the same at 
about 6 mgal RMS for both areas. For comparison, contour maps of the Black and Caspian Sea 
gravity anomalies derived from Rapp’s 180 x 180 geopotential reference model are shown in 
Figures 22 and 23. respectively. Note the close agreement between reference and estimated 
gravity anomalies in the Black Sea. in contrast to the Caspian Sea results. 

Initial estimates of 0.5- x 0.5- gravity anomalies in the Caspian Sea from point geoid data 
on a 0.25 grid were disappointing. The covariance functions were those used in the previous 
study, developed by Rapp from his 180 x 180 gravity reference held. Not only did the estimated 
gravity field look very different from the reference field, despite the close resemblance between 
the geoid height data and the reference model geoid. but. subsequent inverse collocation from 
gravity anomalies to geoid heights did not reproduce the gridded geoid heights. The previous 

study s estimates of gravity anomalies also appear to be very different from those of the 
reference field model. 

A review of the literature for covariance functions indicated that the correlation distance 
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Figure 17. Normalized Rapp's 180 x 180 global covariance functions. 
The normalization coefficients for: 

1) undulation-undulation (solid line) Is 1.13 m 2 , 

2) undulation-gravity (dash line) is 7.26 m-mgal, and 

3) gravity-gravity (dotted line) is 98.36 (mgal) 2 . 
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Rapp’s 180 x 180 global covariance functions. 
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Figure 19. A contour map of the square root of the variances (mgal) 
of the estimated Black Sea gravity anomalies 
based on Rapp's 180 x 180 global covariance functions. 
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Rapp s 180 x 180 global covariance functions. 
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Figure 21. A contour map of the square root of the variances (mgal) 
of the estimated Caspian Sea gravity anomalies 
based on Rapp’s 180 x 180 global covariance functions. 
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Figure 22. A contour map of Rapp s 180 x 180 reference gravity anomalies (mgal) 

of the Black Sea. 
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Figure 23. A contour map of Rapp’s 180 x 180 reference gravity anomalies (mgal) 

of the Caspian Sea* 
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of Rapp’s geoid undulation covariance function is much too long at 3 arc degrees. Models by 
Knudsen ( 1987) and by Jordan Indicated a correlation distance for the geoid beyond degree 180 
of 0.33 and 0.45 arc degrees, respectively. Correlation distance parameters for the gravity 
covariance function and the cross-covariance function also appeared too long, but not by so 
great a factor. As a trial, Jordan's self-consistent set of covariance functions was used, setting 
the geoid function correlation length to 0.5 arc-degree, with the result that the estimated 
gravity anomalies have a closer resemblance to the reference gravity model, and the inverse 
transformation yielded the original geoid. 

Two parameters, C Q and a correlation length, L, are required to calculate Jordan's 
covariance function. Both C Q and L for the undulation-undulation case can be obtained from 
the local empirical covariance function derived from the weighted-average results. We have, 
however, no a-prlort observed gravity anomaly data. A self-consistent approach to estimate 
the Jordan geoid-gravity and gravity-gravity covariance functions would be to empirically 
compute a covariance function from gravity anomalies predicted using Rapp's covariance 
functions. Local empirical covariance functions for the gravity-gravity case can be constructed 
to yield the required Cq and correlation length, L, for Jordan's covariance functions. Program 
JORDAN (Au et ai., 1989) was used to determine consistent sets of Jordan's covariance 
functions. Plots of Jordan's covariance function for the Black Sea are shown in Figure 24. 
Contour maps of the predicted gravity anomalies and square root of the variance for the Black 
Sea are shown in Figures 25 and 26. Note that the predicted gravity error using the Jordan 
covariance function is slightly, but systematically smaller than that of the Rapp covariance 
function in Figure 19. The Jordan covariance functions for the Caspian Sea are shown in Figure 
27; whereas results for the Caspian Sea are shown in Figures 28 and 29. In this case, Jordan's 
covariance function yields a predicted gravity error that is systematically and significantly 
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Figure 24. Normalized Jordan’s theoretical covariance functions of the Black Sea. 
The normalization coefficients for: 

1) undulation-undulation (solid line) is 1.61 m 2 , 

2) undulation-gravity (dash line) is 1 1.61 m-mgal, and 

3) gravity-gravity (dotted line) is 125.44 (mgai) 2 . 
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Figure 26. A contour map of the square root of the variances (mgal) 
of the estimated gravity anomalies of the Black Sea 
based on Jordan’s theoretical covariance functions. 
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Figure 27. Normalized Jordan s theoretical covariance functions of the Caspian Sea. 
The normalization coefficients for: 


1) undulation-undulation (solid line) is 8.47 m 1 2 3 . 

2) undulation-gravity (dash line) is 43.81 m-mgal. and 

3) gravity-gravity (dotted line) is 340 03 (mgal) 2 . 
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Figure 28. A contour map of estimated gravity anomalies (mgal) of the Caspian Sea 
based on Jordan’s theoretical covariance functions. 
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Figure 29. A contour map of the square root of the variances (mgal) 
of the estimated gravity anomalies of the Caspian Sea 
based on Jordan's theoretical covariance functions. 
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larger at 7+ mgal. than that estimated by Rapp’s covariance function (Figure 2 1). These 
differences are apparently due to differences in un-normallzed amplitudes in the respective 
covariance functions. 

The geold-to-gravity transformation is also performed with local empirical covariance 
functions. The empirical covariance functions for the Black Sea are shown in Figure 15. 
Contour maps of gravity anomalies and the square root of the variance for the Black Sea are 
shown In Figures 30 and 31. The corresponding covariance functions for the Caspian Sea are 
shown in Figure 16. Corresponding results for the Caspian Sea are shown in Figures 32 and 
33. The predicted gravity errors are significantly larger at 7+ and 10+ mgal in the Black and 
Caspian Seas, respectively, than those predicted by the Rapp or Jordan covariance functions. 
As before, these increased errors are apparently due to subde differences in the amplitudes of 
the covariance functions. 


51 



Figure 30. A ecu, our map grav lty aucuaMea , mga „ „ *. B,ac k Sea baaeO „„ 

local residual empirical covariance functions. 


52 




BOGHhH^P 



Figure 31. A contour map of the square root of the variances (mgal) 
of the estimated gravity anomalies of the Black Sea 
based on local residual empirical covariance functions. 
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Figure 32. A contour map of estimated gravity ano 
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Figure 33. A contour map of the square root of the variances (mgal) 
of the estimated gravity anomalies of the Caspian Sea 
based on local residual empirical covariance functions. 
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V. DISCUSSION 


A literature search on collocation techniques and local covariance functions for the geoid 
and gravity anomalies reveals that local covariance functions vary markedly from one area to 
another. It should not be a surprise that a single covariance function does not perform equally 
well for both the Black and Caspian Seas. However, one cannot define a local covariance 
function for a given area without detailed knowledge of the gravity and geoid heights in the 
area. This classic "chicken-and-egg" dilemma requires further study. The empirical covariance 
functions developed for the Black and Caspian Seas (Figures 15 and 16) yield qualitatively good 
results in that the resulting predicted gravity anomaly map (Figure 30) closely resembles the 
reference model (Figures 22 and 23) but shows fewer short-wavelength features than the maps 
predicted using the Jordan or Rapp covariance functions (Figures 25 or 18, and 28 or 20). 
Comparison of the shapes of these local covariance functions with those of Rapp (Figure 17) 
or Jordan (Figures 24 and 27) reveal some significant differences. In general, the normalized 
undulation-gravity covariance function exhibits significantly longer correlation distance (point 
on the abscissa where the function drops to 0.5) than the undulation-gravity function, which in 
turn has a longer correlation distance than the gravity-gravity function. This is to be expected 
due to the integral relationship of geoid undulations to gravity. This relationship holds for 
the local empirical covariance functions for the Black Sea (Figure 15). but not for the Caspian 
Sea (Figure 16). The reason for these anomalous properties is not known exactly, but it is 
suspected that the quality of the reference model in the vicinity of the Caspian Sea is of 
critical importance. We note that the empirical covariances for the Black Sea behave normally, 
and there is close agreement between the observed geoid and the model geoid over the Black 
Sea, in contrast to the striking differences between the observed and model geoids over the 
Caspian Sea. Over the Black Sea, accurate, short-wavelength geopotential information has been 
incorporated in Rapp s 180 x 180 geopotential reference model. The altimeter measurements, 
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represented by the gridded geoid heights, add little new information. This is not the case, 
however, for the Caspian Sea. It can be concluded that accurate high-frequency information 
is missing from Rapp s 180 x 180 geopotential model over the Caspian Sea. The quality of the 
transformation from geoid undulations to gravity anomalies apparently varies depending on the 
frequency content of the reference geopotential model, calling into question the robustness of 
the transformation. A consistency test was designed to study the sensitivity of the geoid-to- 
gravity transformation on the information content (or quality) of the reference model. 

A degraded reference model over the Black Sea is obtained by including only long- 
wavelength (36 x 36) terms of Rapp’s 180 x 180 model. Contour maps of the geoid undulations 
and gravity anomalies of the degraded reference model are shown in Figures 34 and 35, 
respectively. The geoid-to-gravity transformation is performed based on the degraded reference 
model and Rapp’s 180 x 180 covariance function. The resultant gravity anomalies are shown in 
Figure 36. It is apparent, comparing with Figure 30, that the quality of the geold-to-gravlty 
transformation is sensitive to the quality of the reference model. A contour map of the 
difference between the estimated gravity anomalies using Rapp’s 180 x 180 reference model and 
36 x 36 reference model is shown in Figure 37. The RMS of the difference is 1 1.35 mgal, with 
several broad areas where the difference exceeds 15 mgal. 

An Iterative transformation to improve the reference model has also been attempted. To 
evaluate this algorithm in a controlled test, Rapp’s 300 x 300 geopotentiai model is adopted as 
the true representation of the geoid and gravity over the Black Sea. and the derived geoid 
surface provides a grid of input data for the geoid-to-gravity transformation and both the 
180 x 180 and 36 x 36 models are used as reference surfaces. The iterative method consists of 
the following steps: 
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Figure 34. A contour map of Rapp’, 36 , 36 reference geotd unduladon, 
(m above mean sea level) of the Black Sea. 
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Figure 35. A contour map of Rapp’s 36 x 36 reference gravity anomalies (mgal) 

of the Black Sea. 
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Figure 36. A eon, our ma p or «!**«, graruy anonyms , mga „ „f , he Black Sra ^ 
Rapp’s 36 x 36 reference geopotential model. 
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Figure 37. A contour map of the difference (mgal) between estimated gravity anomalies 
using Rapp's 180 x 180 and 36 x 36 reference models. 


61 



a) Use Rapp’s covariance function to transform the geoid data to gravity anomalies. 

b) Use the same covariance function to transform the gravity anomalies back to geoid 
heights. 

c) Calculate the RMS difference between the transformed result and the reference model 
for both the geoid and gravity anomalies. 

d) Use the transformed gravity anomalies and transformed geoid results to form a new 
set of local residual empirical covariance functions 

e) Repeat steps a-d using the newly constructed local residual empirical covariance 
functions for transformation and using the transformed results as the new starting 
reference model, until the RMS difference satisfies a convergence criterion. 

It is observed that using the iterative process for both the 180 x 180 and 36 x 36 reference 
models, a major correction to the reference models occurs during the first and second Iterations. 
The Iterative process converges in less than five iterations for reasonable Integration cap radius, 
such as one degree. The transformed gravity anomalies based on 180 x 180 and 36 x 36 
Reference models are shown in Figures 38 and 39. respectively. Contour maps of the "true" 
geoid undulations and gravity anomalies according to Rapp’s 300 x 300 model are shown in 
Figures 40 and 41, respectively. The iterative method does not seem to materially improve the 
transformation when a poor (36 x 36 compared to 180 x 180) reference surface is used. 

The sensitivity of the geoid-to-gravity transformation to different covariance functions and 
information content of the reference models is quantified by determining the RMS of the 
difference between the "true" gravity anomalies and the estimated ones, as shown in Table 1. 

The RMS values represent the error of commission in the geoid-to-gravity transformation. To 
make the comparison fair, a common a-prtorl stabilizing variance of |25 cm| 2 was added for all 
three covariance functions when the 180 x 180 reference field was used, and (60 cm| 2 was added 
whenever the 36 x 36 reference field was used. The iterative algorithm based on empirical 
covariance functions, compared to the single-pass transformation based on Rapp’s and Jordan s 
covariance functions, generally yields the best recovered gravity anomalies when the information 
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Figure 38. A contour map of the estimated gravity anomalies (mgal) of the Black Sea 
based on Rapp’s 180 x 180 reference geopotential model 
and the seif-consistent iterative approach. 

Rapp’s 300 x 300 model geoid undulations are used as input data. 

This is a test of the iterative approach. 
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Figure 39. A contour map of the estimated gravity anomalies (mgal) of the Black Sea 
based on Rapp's 36 x 36 reference geopotential model 
and the self-consistent iterative approach. 

Rapp's 300 x 300 model geoid undulations are used as input data. 

This is a test of the iterative approach. 
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Figure 40. A contour map of Rapp's 300 x 300 reference geoid undulations 
(m above mean sea level) of the Black Sea. 
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Figure 4 1. A contour map of Rapp’s 300 x 300 reference gravity anomaiies (mgal) 

of the Black Sea. 
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spectra limit of the reference model Is commensurate with the cap size ( 180 x 180 model and 
1° cap size). 


transformation 

covariance 

functions 

180 x 180 
reference model 
(mgal RMS) 

36x36 

reference model 
(mgal RMS) 

Rapp's 

4.47 

14.39 

Jordan’s 

4.34 

11.12 

Iterative 

3.19 

11.79 


Table 1 Error of commission in the geoid-to-gravity transformation as a function 
of different transformation algorithms using 1° cap size for the integration 
region. Rapp's 180 x 180 and 36 x 36 models are used as reference 
surfaces. 

When the stabilizing a-priort error variances are removed or replaced by lower values, 
while still maintaining solution stability, the RMS error of commission is generally reduced, 
along with the formal prediction error. For example, in the case of Jordan’s covariance 
function, when the a-priort variance is reduced to ( 10 cm) 2 , the corresponding values in Table 
1 become 2.97 mgal and 6.96 mgal respectively for 180 x 180 and 36 x 36 reference models. 

In fact, these levels of stabilizing noise seem to be optimal for both the Jordan and Rapp 
covariance functions, because lower and higher values of a-prtorl noise result in larger RMS 
errors of commission. The error of commission for Rapp’s covariance function Is universally a 
few percent higher than that for Jordan’s covariance function. On the other hand, (60 cmp is 
optimal (and necessary) for the empirical local covariance function. Based on these results, 
we conclude that Jordan's covariance function is best for gravity prediction in the Black Sea 
region, and that the optimal level of a-prtorl noise is about 10 cm. Furthermore, we conclude 
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that the error of commission of the least-squares collocation technique for gravity prediction 
is highly dependent on the quality of the reference model, ranging from 3 to 7 mgal. 

It should be appreciated that collocation is a statistical method that relies on the 
transformation covariance functions to provide the physics of the figure of the Earth and its 
gravity field. The iterative algorithm is a hybrid of a perturbation on the reference surface 
and an information shaping / filtering process. The initial shaping filter, the covariance 
function, should conform with the information content of the initial reference surface. That 
is, for a reference model whose short wavelength cutofT is at 1°, the covariance function should 
represent Information of wavelengths shorter than 1° and the integration cap size should have 
a commensurate size. It is speculated that when the 36 x 36 model is used as a reference 
surface, long wavelength correction to the updated reference surface is limited to wavelengths 
less than 1° because of the chosen integration cap size of 1°. Wavelength components longer 
than 1° must be corrected in order to improve the 36 x 36 reference model. The integration 
cap size, therefore, should be relatively large if a long wavelength reference model is to be 
used. A large integration cap size, unfortunately, will result in forbiddingly high computing 
costs unless the data grid for the initial iterative steps is decimated. The cap size can be 
gradually reduced, as the data density is gradually increased, in subsequent iterations. However, 
an algorithm developed to maintain constant density of data in each iterative step, performed 
worse than the conventional algorithm. 
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